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We perform numerical simulations to examine particle diffusion at steady shear in a model gran- 
ular material in two dimensions at the jamming density and zero temperature. We confirm findings 
by others that the diffusion constant depends on shear rate as D ~ 7 9D with qu < 1, and set out 
to determine a relation between qo and other exponents that characterize the jamming transition. 
We then examine the the velocity auto-correlation function, note that it is governed by two pro- 
cesses with different time scales, and identify a new fundamental exponent, A, that characterizes an 
algebraic decay of correlations with time. 

PACS numbers: 45.70.-n, 64.60.-i 



As the volume fraction increases in zero-temperature 
collections of spherical particles with repulsive contact 
interaction, there is a transition from a liquid to an amor- 
phous solid state — the jamming transition. It has been 
suggested that this transition is a critical phenomenon 
with universal critical exponents [1] and the properties 
of this transition continues to be a very active field of 
research. Simulations at steady shearing have provided 
strong evidence that the behavior at the jamming den- 
sity actually is a critical phenomenon [2, 3], but questions 
still remain to what extent results and ideas from ordi- 
nary critical phenomena may be taken over to the study 
of jamming as well as the fundamental reason for the 
observered critical behavior. 

In critical phenomena the behavior is governed by a di- 
verging length scale and one expects that this should also 
be reflected in the time dependence of various quantities. 
One way to probe the time dependence is to measure 
the particle displacements and thereby the diffusion con- 
stant. Experiments suggest that the diffusion depends 
algebraic on the shear rate, D ~ j qD , with qo < 1 [4, 5]. 
Since this appears to be one more critical exponent, and 
one usually expects relations between different critical ex- 
ponents, the existence of such a relation between qu and 
other exponents that characterize the jamming transition 
is an interesting question. 

In this Letter we examine the velocity auto-correlation 
function in an attempt to understand the behavior of 
the diffusion constant. A careful study of this function 
at very low shear rates reveals that it has both an al- 
gebraic decay and an exponential cutoff. It futhermore 
turns out that these two processes are governed by two 
different time scales, with the exponential decay being re- 
lated to the externally applied time scale ~ 7 _1 whereas 
the remaining part — which we identify with an internal 
relaxation — is governed by a time scale ~ er -1 . 

Following O'Hern et al. [6] we simulate frictionless soft 
disks in two dimensions using a Bi-dispersive mixture 
with equal numbers of disks with two different radii of 
ratio 1.4. Length is measured in units of the small par- 
ticles (d s — 1). With rij for the distance between the 



centers of two particles and dij the sum of their radii, 
the interaction between overlapping particles is 

T/fV \ — f 2(1 ~ r ij / dij) 2 1 r ij < dij, 
''' \<). Tij > dij. 

We use Lees-Edwards boundary conditions [7] to intro- 
duce a time-dependent shear strain 7 = tj. With peri- 
odic boundary conditions on the coordinates Xi and yi 
in an L x L system, the position of particle i in a box 
with strain 7 is defined as = (xi + yyi, yi). We simu- 
late overdamped dynamics at zero temperature with the 
equation of motion [8] , 

drj _ ^ dVjrij) 

The unit of time is To = d s /Ce. We take e = 1 and 
C = 1. We integrate the equations of motion with the 
Heuns method, using a time step At = 0.2to. As this 
must be considered rather large, we have checked care- 
fully that simulations with half that time step gives the 
same results to a very high accuracy. The possibility to 
use such large time steps is linked to the simple dynamics, 
zero temperature, and our low shear rates. 

We study a system with many particles, N = 65536, at 
4>j, since the correlation length in the system should only 
depend on the finite shear rate and one therefore expects 
a simpler behavior. We have checked that our results are 
not affected by finite size effects. The behavior of the 
shear stress at three densities at and around (j> — 0.8428 
is shown in Fig. 1. At <f> = 0.8428, which is our candidate 
for (f)j, the shear stress is algebraic in the shear rate, a ~ 
j q " with q„ — 0.386, whereas the data away from <f>j have 
clear curvatures. In the notation of Ref. [2], q a = A/(/? + 
A). We remark that the fit is not entirely perfect, in 
spite of the nice algebraic behavior in Fig. 1. This is the 
reason why the present estimate (j>j = 0.8428 is somewhat 
higher than the estimate in Ref. [2] . Our new data (which 
extends down to lower shear rates) also show that a high- 
precision determination of 4>j and the related exponents 
is a difficult task. This is due to some corrections to the 
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FIG. 1: Behavior at <j>j. At cj> = 0.8428 which is a good 
candidate for <j)j the shear stress depends algebraically on j 
to a very good approximation, a ~ A / q " , with g CT = 0.386. 
Data at higher (squares) and lower (crosses) densities show 
clear curvatures. 



expected scaling behavior, as will be discussed elsewhere. 
For the purpose of the present Letter the approximate 
value 4>j « 0.8428 is, however, entirely sufficient. 

We determine the diffusion constant from the trans- 
verse displacements, i.e. the displacements in the y di- 
rection, and the velocity auto-correlation function from 
the y component of the velocity, 

9v (t) = (Vy(t')Vy(t' +<)), 

where the average is over all particles and a large number 
of initial times, t' . Here and in the following, t is the 
difference between two absolute times. The velocity auto- 
correlation function has been examined before [9] , but the 
present data with higher precision at lower shear rates 
makes it possible to do a more thorough analysis of its 
properties. The relation to the diffusion constant is given 
by the fundamental relation 



D = 



dtg v (t) =g v (0) / dtG v (t) 
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FIG. 2: Particle displacements for 7 = 10~ 6 . Panel (a) shows 
the crossing over of (Ay 2 ) from ballistic behavior at short 
times to diffusion at large times. The slopes of the solid lines 
are 2 and 1, respectively. Panel (b) shows the probability dis- 
tribution function of particle displacements, normalized by 
the width of the respective distributions. Note the exponen- 
tial shape at small strains (= short times) that crosses over 
to a Gaussian distribution at 7 ~ 0.2. Panel (c) shows the 
determination of D from the large 7 part of the same data. 
The solid line is from fitting (Ay 2 ) to Eq. (2); the dashed 
line corresponds to D. 



where we introduce the normalized G v (t) = g v (t)/ g v (0). 
It is convenient to write the expression in terms of G v (t) 
both since it is the quantity that will be examined be- 
low and since the prefactor, g v (0), has a known behav- 
ior, g v (0) = Vy ~ cry ~ ji+i" 5 which follows from 
7V(v 2 )/C = L 2 cr 7 [9]. 

Some quantities related to the particle displacements 
are shown in Figs. 2. To make it easier to interpret the 
figures these quantities are plotted against 7 (the strain 
increment), though we discuss the behavior in terms of 
t. Panel (a) which shows (Aj/ 2 ) against 7, illustrates the 
crossover from ballistic motion at short times to diffu- 
sion, (Ay 2 ) ~ t. The probability distribution function 
(PDF) of Ay (normalized by the width of the distribu- 
tion), for several different strain increments, is shown in 



panel (b). The PDF crosses over from exponential be- 
havior at short times (small 7) to a Gaussian at longer 
times, as found by others [5, 10]. Our determination of 
the diffusion constant is illustrated in Fig. 2(c). As the 
figure shows it is difficult to determine D from the long 
time limit of (Ay 2 ) jt since this quantity approaches the 
constant value = D very slowly. The reason for this is a 
remainder of the short time behavior. For t > t$, where 
to is the range of the velocity correlations (such that 
G v {t) may be neglected for t > t ; we choose 70 = 0.5, 
to = 7o H) ■, h is easy to show that the expression for the 
mean square distance is 

(Ay 2 (t)) = f dt' [ dt"g v (t' - t") =Dt-d (2) 
Jo Jo 
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FIG. 3: Diffusion constant versus shear rate. The open circles 
are the diffusion constant at <j> = 0.8428 versus 7. The data 



is well fitted to an algebraic relation D 
0.785(5). 



with qo 



with D from Eq. (1) and d = /„" dt' f t ° dt"g v {t"). The 
solid line in Fig. 2(b) is from a fit to Eq. (2) with data 
from the interval 7o = 0.5<7<2. The dashed line is 
the estimated value of D. 

Figure 3 shows diffusion constant versus shear rate, 
determined with the same kind of fits. The behavior is 
D ~ 7« D , with q D = 0.785(5). This implies that the 
distance moved per unit strain decreases with increasing 
shear rate [11]. The corresponding exponents from ex- 
periments are qD = 0.80 ± 0.01 from three dimensional 
colloids [4] and qo = 0.66 ± 0.05 from bubble rafts [5]. 
We note that the experiments on the colloids were per- 
formed at a density close to the jamming density whereas 
the bubble raft was studied well above 4>J- This is a pos- 
sible reason why the value of qn in the colloids agrees 
well with our value obtained at <pj. 

To examine this behavior we turn to the velocity auto- 
correlation function which is shown in Fig. 4(a) for a 
range of shear rates. The same data is shown also in panel 
(b), but now plotted against y(t) = tj with a linear scale 
on the x axis. From this figure it seems that logG.„ at 
large 7 behaves linearly with similar slopes for different 
7, which suggests an exponential decay, ~ e ~ 7 ^^ 71 . Wc 
take this to suggest that G„(i,7) may be written 

G v (t,'j) = G i » t (t, 7 ) e-*^\ 



which means that G v is a product of an exponential decay 
governed by the externally imposed time scale t\ =71/7 
and a function GJ, nt , which captures the internal relax- 
ational dynamics. 

In the attempt to make sense of this data we found 
that a certain choice of 71 in Eq. (3) leads to a great sim- 
plification. It turns out that GJf t (i,7) is then a function 
of the combination ty K , 

Gf(t, 7 ) = G(n K ), 

and, furthermore, that this function at large values of 
its argument behaves algebraically, G(x) ~ x~ x , see 
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FIG. 4: The velocity auto-correlation function. Panel (a) 
shows the normalized auto-correlation function versus time 
and panel (b) is the same data against strain f(t). The key 
to our further analysis is the rectilinear behavior at large 7 
in panel (b) which suggests an exponential decay at large 7, 
~ e 'w/t 1 . The line corresponds to 71 = 0.06. Panel (c) 
is the same data but now compensated for such an exponen- 
tial decay with 71 = 0.088 (from fitting to Eq. (4)). Panel 
(d) shows the collapse, now with to as the scaling variable. 
The scaling function obeys G(x) ~ x~ x for large values of its 
argument. 



(3) Fig. 4(c). This defines a new fundamental exponent that 



characterizes the relaxation dynamics, 



r 



though- 



for the accessible shear rates — it is to some extent masked 
by the exponential decay. To get unbiased values for 
these parameters, we fitted all data in the range 0.002 < 
G„(t,7) < 0.3 to 



G,(t, 7 ) 
A, 



a (tr) 



(4) 



with A, 71, A, and n as free parameters. We find 
71 = 0.088, A = 0.78, and k = 0.384. The solid lines 
in Fig. 4(c) show ^4(^7 K )~ A for different 7. 

In this expression, assumes the role of a charac- 
teristic time for the internal relaxation. The observation 
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that Ca has the dimension of inverse time together with 
the good numerical agreement between k — 0.384 and 
q a = 0.386 (recall that a ~ y°) suggests that y can be 
substituted with a such that the scaling function may be 
written G(ta). Panel (d) shows the collapse when plot- 
ting against the scaling variable ta. Note that 71 is the 
only adjustable parameter in this plot. 

We now like to determine the relation between qjj and 
the exponents q„ and A that characterize the scaling of 
£?.„(£, 7). In the first approximation we neglect the satu- 
ration of G v at small t, and, in effect, assume that Eq. (4) 
holds down to t = 0. This gives 



D 



A-l 



dx x 



7 



which leads to q^ = q a + \ — n\ = A+(l — \)q a = 0.865. 
This is the expected behavior as 7 — ► 0, but since it is 
derived from a simplified G v (t,j) and the result is well 
above qo — 0.785 from Fig. 3, we next try to take the 
saturation of G v (t, 7) at small t into account and write 




FIG. 5: Behavior of the diffusion constant. In the limit of low 

shear rate we expect the behavior which is given by the solid 

(i) 

line, D = A\^ q o . The open circles include the corrections to 
this behavior as given by Eq. (5). Note the good agreement 
with the measured D from Fig. 3. We conclude that qn ~ 
0.785 from the solid symbols (cf. Fig. 3), is only an effective 
exponent that describes the behavior in a limited range of 
shear rates. 
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Assuming that e~*^ 71 rj I at t — £o/7 K (which holds 
to a good approximation up to our highest shear rate, 
7 = 10~ 5 ) the function is continuous at £0 if A = £q . We 
then get 



D = A x yo - A 27 , 



with A\ = 0.235 and A 2 = 0.480. As shown in Fig. 5 
this expression (open circles) approaches D ~ y° (solid 
line) at small 7 whereas there is an appreciable difference 
with a smaller slope at larger 7. The diffusion constant 
from Fig. 3 is shown as solid dots. Note that it agrees 
well with the open circles from Eq. (5). It therefore seems 
that the measured exponent qjj = 0.785 is not the true 
asymptotic behavior. 

A central conclusion from our analysis is that the full 
G v (t,j) is approaching an algebraic behavior ~ t~ x as 
7^0. We now speculate that the algebraic behav- 
ior is related to the finding from quasistatic simulations 
that individual plastic events often are avalanches of el- 
ementary flips [12-15]. The reason for making this con- 
nection is ideas from self-organized criticality — with the 
paradigmatic sandpile model — that a driven system can 
automatically adjust itself such that there are avalanches 
on all length and time scales, which would be seen 
through power laws. With a sufficiently low shear rate 
(say 7 = 10 -10 or 10~ 9 ) there would be time for the 
avalanches to occur one at a time and evolve according 
to their own dynamics. At higher shear rates other ef- 
fects appear that kill off the avalanches. One possible 
mechanism is that a new avalanche interferes with an 



existing one and thereby destroys its internal dynamics. 
Another possibility is that it is simply the shearing of the 
simulation box that destroys the correlations. 

To conclude, we have found that the velocity auto- 
correlation function is governed by two different time 
scales. With t\ — 71/7 from the externally applied 
shear rate and ii„t = ~ r y~ qa for the internal relax- 
ation, the velocity auto-correlation function is G v (t, 7) = 
(5) G(i/ii n t)e _t / tl , where G(x) ~ x~ x for large x and A is a 
new fundamental exponent. This also leads to the desired 
expression for qjj in terms of two fundamental exponents, 
qo = A + (1 — X)q a . We speculate that this algebraic de- 
cay is related to avalanches of elementary flips, and could 
be a manifestation of self-organized criticality. 
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